#include <string.h>
#include <cmath>
#include <iostream>
#include <sstream>
#include <vector>
#include <stack>
#include <list>

#include "matrix.h"

#define USING_LIST_AS_STACK

using namespace std;

jieningMatrix::jieningMatrix(unsigned int xlen, unsigned int ylen, double initvalue)
: m_xlen(1), m_ylen(1)
{
	if(0 < xlen || xlen<= 10000)
	{
		this->m_xlen = xlen;
	}

	if(0 < ylen || ylen <= 10000)
	{
		this->m_ylen = ylen;
	}

	this->m_values = new double[this->m_xlen*this->m_ylen];
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] = initvalue;
	}
}

jieningMatrix::jieningMatrix(unsigned int xlen, unsigned int ylen, double initvalues[])
: m_xlen(1), m_ylen(1)
{
	if(0 < xlen || xlen<= 10000)
	{
		this->m_xlen = xlen;
	}

	if(0 < ylen || ylen <= 10000)
	{
		this->m_ylen = ylen;
	}

	this->m_values = new double[this->m_xlen*this->m_ylen];
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] = initvalues[i];
	}
}

jieningMatrix::jieningMatrix(const jieningMatrix& obj)
{
	//cout<<"in jieningMatrix::jieningMatrix(const jieningMatrix& obj)"<<endl;
	this->m_xlen = obj.m_xlen;
	this->m_ylen = obj.m_ylen;
	this->m_values = new double[this->m_xlen*this->m_ylen];
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] = obj.m_values[i];
	}
}

//赋值符重载
jieningMatrix& jieningMatrix::operator =(const jieningMatrix& obj)
{
	if(this == &obj)
	{
		return *this;
	}

	delete this->m_values;

	this->m_xlen = obj.m_xlen;
	this->m_ylen = obj.m_ylen;
	this->m_values = new double[this->m_xlen*this->m_ylen];
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] = obj.m_values[i];
	}
	return *this;
}

//加号重载
jieningMatrix jieningMatrix::operator +(const jieningMatrix& obj)
{
	if(this->m_xlen != obj.m_xlen || this->m_ylen != obj.m_ylen)
	{
		return *this;
	}

	jieningMatrix rtnobj(this->m_xlen, this->m_ylen);
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		rtnobj.m_values[i] = this->m_values[i] + obj.m_values[i];
	}
	return rtnobj;
}

jieningMatrix& jieningMatrix::operator +=(const jieningMatrix& obj)//加号重载
{
	if(this->m_xlen != obj.m_xlen || this->m_ylen != obj.m_ylen)
	{
		return *this;
	}

	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] += obj.m_values[i];
	}
	return *this;
}

jieningMatrix jieningMatrix::operator +(double val)
{
	jieningMatrix rtnobj(this->m_xlen, this->m_ylen);
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		rtnobj.m_values[i] = this->m_values[i] + val;
	}
	return rtnobj;
}

jieningMatrix& jieningMatrix::operator +=(double val)
{
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] += val;
	}
	return *this;
}

//减号重载
jieningMatrix jieningMatrix::operator -(const jieningMatrix& obj)
{
	if(this->m_xlen != obj.m_xlen || this->m_ylen != obj.m_ylen)
	{
		return *this;
	}

	jieningMatrix rtnobj(this->m_xlen, this->m_ylen);
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		rtnobj.m_values[i] = this->m_values[i] - obj.m_values[i];
	}
	return rtnobj;
}

jieningMatrix& jieningMatrix::operator -=(const jieningMatrix& obj) //减号重载
{
	if(this->m_xlen != obj.m_xlen || this->m_ylen != obj.m_ylen)
	{
		return *this;
	}

	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] -= obj.m_values[i];
	}
	return *this;
}

jieningMatrix jieningMatrix::operator -(double val)
{
	jieningMatrix rtnobj = *this;
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		rtnobj.m_values[i] = this->m_values[i] - val;
	}
	return rtnobj;
}

jieningMatrix& jieningMatrix::operator -=(double val)
{
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] -= val;
	}
	return *this;
}
//乘号重载
jieningMatrix jieningMatrix::operator *(const jieningMatrix& obj)
{
	if(this->m_xlen != obj.m_ylen)
	{
		return *this;
	}
	jieningMatrix rtnobj(obj.m_xlen, this->m_ylen);
	unsigned int loop=0;
	for(unsigned int y=0; y<this->m_ylen; ++y)
	{
		for(unsigned int x=0; x<obj.m_xlen; ++x)
		{
			rtnobj.m_values[loop] = 0.0;
			double tmpv;
			for(unsigned int i=0; i<this->m_xlen;++i)
			{
				tmpv = this->getValue(i, y);
				tmpv *= ((jieningMatrix&)obj).getValue(x, i);
				rtnobj.m_values[loop] += tmpv;
			}
			++loop;
		}
	}
	return rtnobj;
}

jieningMatrix& jieningMatrix::operator *=(const jieningMatrix& obj) //乘号重载
{
	if(this->m_xlen != obj.m_ylen)
	{
		return *this;
	}
	jieningMatrix rtnobj(obj.m_xlen, this->m_ylen);
	unsigned int loop=0;
	for(unsigned int y=0; y<this->m_ylen; ++y)
	{
		for(unsigned int x=0; x<obj.m_xlen; ++x)
		{
			rtnobj.m_values[loop] = 0.0;
			double tmpv;
			for(unsigned int i=0; i<this->m_xlen;++i)
			{
				tmpv = this->getValue(i, y);
				tmpv *= ((jieningMatrix&)obj).getValue(x, i);
				rtnobj.m_values[loop] += tmpv;
			}
			++loop;
		}
	}

	this->m_xlen = rtnobj.m_xlen;
	this->m_ylen = rtnobj.m_ylen;
	delete this->m_values;
	this->m_values = new double[this->m_xlen*this->m_ylen];
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] -= rtnobj.m_values[i];
	}

	return *this;
}

jieningMatrix jieningMatrix::operator *(double val)
{
	jieningMatrix rtnobj = *this;
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		rtnobj.m_values[i] = this->m_values[i] * val;
	}
	return rtnobj;
}

jieningMatrix& jieningMatrix::operator *=(double val)
{
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] *= val;
	}
	return *this;
}

double jieningMatrix::getValue(unsigned int xpos, unsigned int ypos)
{
	if(xpos>=this->m_xlen)
	{
		return 0.0;
	}
	if(ypos>=this->m_ylen)
	{
		return 0.0;
	}
	return this->m_values[ypos*this->m_xlen + xpos];
}

string jieningMatrix::toString() const
{
	ostringstream strbuf;
	strbuf<<endl;
	unsigned int loop=0;
	for (unsigned int y = 0; y < this->m_ylen; ++y)
	{
		strbuf<<y<< " --> | ";
		for (unsigned int x = 0; x <  this->m_xlen; ++x)
		{
			strbuf <<x<<"->"<<this->m_values[loop]<<"\t";
			++loop;
		}
		strbuf<<" |" << endl;
	}
	strbuf<<"matrix["<<this->m_xlen<<" * "<<this->m_ylen<<"]"<<endl;
	//strbuf<<endl;
	return string(strbuf.str());
}

void jieningMatrix::doPrint() const
{
	cout<<this->toString();
	return;
	cout<<endl;
	unsigned int loop=0;
	for (unsigned int y = 0; y < this->m_ylen; ++y)
	{
		cout <<y<<"| ";
		for (unsigned int x = 0; x <  this->m_xlen; ++x)
		{
			cout <<this->m_values[loop]<<"\t";
			++loop;
		}
		cout<<" |" << endl;
	}
	cout<<"matrix["<<this->m_xlen<<" * "<<this->m_ylen<<"]"<<endl;
}

// 去掉第y（0 base）行第x（0 base）列后子矩阵
// 如果 xdrop = -1表示返回去掉第y（0 base）行后的子矩阵
// 如果 ydrop = -1表示返回去掉第x（0 base）列后的子矩阵
jieningMatrix jieningMatrix::getSubMatrix(int xdrop, int ydrop)
{
	if((-1 != xdrop && 0 == (this->m_xlen-1)) || 
		(-1 != ydrop && 0 ==  (this->m_ylen-1)))
	{
		return *this;
	}
	if(xdrop >= (int)this->m_xlen || ydrop >= (int)this->m_ylen)
	{
		return *this;
	}
	
	unsigned int rtnmxlen = (-1 == xdrop) ? this->m_xlen : (this->m_xlen-1);
	unsigned int rtnmylen = (-1 == ydrop) ? this->m_ylen : (this->m_ylen-1);

	jieningMatrix rtnobj(rtnmxlen, rtnmylen);
	unsigned int loopr=0;
	unsigned int loopw=0;
	for (int y = 0; y < (int)this->m_ylen; ++y)
	{
		for (int x = 0; x < (int)this->m_xlen; ++x)
		{
			if(x != xdrop && y != ydrop)
			{
				rtnobj.m_values[loopw] = this->m_values[loopr];
				++loopw;
			}
			++loopr;
		}
	}
	return rtnobj;
}

unsigned int jieningMatrix::getXlen()
{
	return this->m_xlen;
}
unsigned int jieningMatrix::getYlen()
{
	return this->m_ylen;
}

bool jieningMatrix::getMData(double* cptoptr, unsigned int cpsize)
{
	if(NULL == cptoptr)
	{
		return false;
	}
	cpsize = (cpsize<=(this->m_xlen*this->m_ylen)) ? cpsize : (this->m_xlen*this->m_ylen);
	memcpy(cptoptr,  this->m_values, cpsize*sizeof(double));
	return true;
}

//交换r1行和r2行
bool jieningMatrix::swapRow(unsigned int r1, unsigned int r2)
{
	if(r1>=this->m_ylen || r2>=this->m_ylen)
	{
		return false;
	}
	double tmpforswap;
	for(unsigned int loop =0; loop<this->m_xlen; ++loop)
	{
		tmpforswap = this->m_values[r1*this->m_xlen + loop];
		this->m_values[r1*this->m_xlen + loop] = this->m_values[r2*this->m_xlen + loop];
		this->m_values[r2*this->m_xlen + loop] = tmpforswap;
	}
	return true;
}

//交换line1列和line2列
bool jieningMatrix::swapLine(unsigned int line1, unsigned int line2)
{
	if(line1>=this->m_xlen || line2>=this->m_xlen)
	{
		return false;
	}
	double tmpforswap;
	for(unsigned int loop =0; loop<this->m_ylen; ++loop)
	{
		tmpforswap = this->m_values[loop*this->m_xlen + line1];
		this->m_values[loop*this->m_xlen + line1] = this->m_values[loop*this->m_xlen + line2];
		this->m_values[loop*this->m_xlen + line2] = tmpforswap;
	}
	return true;
}

// 高斯消元法，解线性方程组
jieningMatrix jieningMatrix::computeFCZ()
{
	jieningMatrix rtnobj(*this);
		
	for(unsigned int row=0;row<rtnobj.m_ylen;++row)
	{
		//cout<<"in computeFCZ row="<<row<<"\t"<<rtnobj->getValue(rtnobj->m_xlen-1, 0)<<endl;
		//rtnobj->doPrint();
		if(row >= rtnobj.m_xlen)
		{
			break;
		}
		//todo 从第row行开始 寻找一行其第row列数字不是0,把这一行交换到第row行
		unsigned int findrow=row;
		bool findnone = true;
		for(;findrow<rtnobj.m_ylen;++findrow)
		{
			if(abs(rtnobj.getValue(row, findrow)) > jieningMatrix::ZEROVAL)
			{
				findnone = false;
			}
		}
		if(findnone)
		{
			// 没有找到这样行，处理下一行
			continue;
		}
		if(findrow != row)
		{		
			// 如果找到的行不是第row行，和row行交换
			rtnobj.swapRow(findrow,row);
		}
		
		// 第row行 第row列 变换为1.0
		double rowrowval = rtnobj.m_values[row*rtnobj.m_xlen + row];
		for(unsigned loop=row+1; loop<rtnobj.m_xlen; ++loop)
		{
			if(abs(rtnobj.m_values[row*rtnobj.m_xlen + loop]) > jieningMatrix::ZEROVAL)
			{
				rtnobj.m_values[row*rtnobj.m_xlen + loop] /= rowrowval;
			}
		}
		rtnobj.m_values[row*rtnobj.m_xlen + row] = 1.0;

		// 消元其余行 第row列
		for(unsigned looprow = 0;looprow<rtnobj.m_ylen; ++looprow)
		{
			if(looprow == row)
			{
				continue;
			}
			else
			{
				double xiaoYuanXiShu = rtnobj.m_values[looprow*rtnobj.m_xlen + row];
				if(abs(xiaoYuanXiShu) > jieningMatrix::ZEROVAL)
				{
					for(unsigned loop=row+1; loop<rtnobj.m_xlen; ++loop)
					{
						double samelineInrow = rtnobj.m_values[row*rtnobj.m_xlen + loop];
						if(abs(samelineInrow) > jieningMatrix::ZEROVAL)
						{
							rtnobj.m_values[looprow*rtnobj.m_xlen + loop] -= xiaoYuanXiShu * samelineInrow;
						}
					}
					rtnobj.m_values[looprow*rtnobj.m_xlen + row]  = 0.0;
				}
			}
		}//end of for()
	}// end of for()
	
	return rtnobj;
}


jieningMatrix::~jieningMatrix()
{
	//cout<<"in jieningMatrix::~jieningMatrix"<<endl;
	delete this->m_values;
	this->m_values = NULL;
}

//============================================================
//============================================================
//============================================================

jieningSquareMatrix::jieningSquareMatrix(unsigned int len, double initvalues[])
{
	if(0 < len || len<= 10000)
	{
		this->m_xlen = this->m_ylen = len;
	}
	else
	{
		this->m_xlen = this->m_ylen = 1;
	}
	this->m_values = new double[this->m_xlen*this->m_ylen];
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] = initvalues[i];
	}
}

jieningSquareMatrix::jieningSquareMatrix(unsigned int len, double initvalue)
{
	if(0 < len || len<= 10000)
	{
		this->m_xlen = this->m_ylen = len;
	}
	else
	{
		this->m_xlen = this->m_ylen = 1;
	}

	this->m_values = new double[this->m_xlen*this->m_ylen];
	unsigned int loop=0;
	for (unsigned int y = 0; y < this->m_ylen; ++y)
	{
		for (unsigned int x = 0; x <  this->m_xlen; ++x)
		{
			if(x == y)
			{
				this->m_values[loop] = initvalue;
			}
			else
			{
				this->m_values[loop] = 0.0;
			}
			++loop;
		}
	}
}

//拷贝构造函数
jieningSquareMatrix::jieningSquareMatrix(const jieningSquareMatrix& obj)
{
	this->m_xlen = this->m_ylen = obj.m_xlen;
	this->m_values = new double[this->m_xlen*this->m_ylen];
	for(unsigned int i=0; i<this->m_xlen*this->m_ylen;++i)
	{
		this->m_values[i] = obj.m_values[i];
	}
}

//拷贝构造函数
jieningSquareMatrix::jieningSquareMatrix(const jieningMatrix& obj)
{
	unsigned int objxlen = ((jieningMatrix&)obj).getXlen();
	unsigned int objylen = ((jieningMatrix&)obj).getYlen();
	this->m_xlen = this->m_ylen = (objxlen < objylen ? objxlen : objylen);
	this->m_values = new double[this->m_xlen*this->m_ylen];

	double* objvalues = new double[objxlen*objylen];
	((jieningMatrix&)obj).getMData(objvalues, objxlen*objylen );

	unsigned int loopr=0;
	unsigned int loopw=0;
	for (unsigned int y = 0; y < objylen; ++y)
	{
		for (unsigned int x = 0; x < objxlen; ++x)
		{
			if(x < this->m_xlen && y <  this->m_ylen)
			{
				this->m_values[loopw] = objvalues[loopr];
				++loopw;
			}
			++loopr;
		}
	}

	delete objvalues;
}

jieningSquareMatrix::~jieningSquareMatrix()
{
	//cout<<"in ~jieningSquareMatrix"<<endl;
}

// 递归计算矩阵
double jieningSquareMatrix::doDetComputR()
{
	if(1 == this->m_xlen)
	{
		return this->m_values[0];
	}
	if(2 == this->m_xlen)
	{
		double tmpv = this->m_values[0]*this->m_values[3] - this->m_values[1]*this->m_values[2];
		return tmpv;
	}

	int fuhao=1;
	double tmpv=0.0;
	for(unsigned int i=0; i< this->m_xlen; ++i)
	{
		if(abs(this->m_values[i]) < jieningMatrix::ZEROVAL)
		{
		}
		else
		{
			jieningSquareMatrix subSM = this->getSubMatrix(i, 0);
			//subSM.doPrint();
			double subdetv = subSM.doDetComputR();
			tmpv += (double)fuhao*this->m_values[i]*subdetv;
			//cout<<"in doDetComput tmpv="<<tmpv<<"\t mv="<<this->m_values[i]<<"\t subdet="<<subdetv<<endl;
		}
		fuhao = -fuhao;
	}
	return tmpv;
}

// 消除递归的计算矩阵
double jieningSquareMatrix::doDetComput()
{
	if(1 == this->m_xlen)
	{
		return this->m_values[0];
	}
	if(2 == this->m_xlen)
	{
		double tmpv = this->m_values[0]*this->m_values[3] - this->m_values[1]*this->m_values[2];
		return tmpv;
	}
	if(3 == this->m_xlen)
	{
		double tmpv = 0.0;
		tmpv += this->m_values[0]*(this->m_values[4]*this->m_values[8] - this->m_values[5]*this->m_values[7]);	
		tmpv -= this->m_values[1]*(this->m_values[3]*this->m_values[8] - this->m_values[5]*this->m_values[6]);	
		tmpv += this->m_values[2]*(this->m_values[3]*this->m_values[7] - this->m_values[4]*this->m_values[6]);	
		return tmpv;
	}


#ifndef USING_LIST_AS_STACK
	//* using stack container
	stack<double> xisuStack;
	stack<jieningSquareMatrix> smStack;
	xisuStack.push(1.0);
	smStack.push(*this);
#else
	//* using list container
	list<double> xisuStack;
	list<jieningSquareMatrix> smStack;
	xisuStack.push_back(1.0);
	smStack.push_back(*this);

#endif

	double tmpv=0.0;
	int loopcount=0;
	while(true)
	{
		++loopcount;
		if(0 == smStack.size())
		{
			break;
		}
		
		cout<<endl<<"in doDetComput loop="<<loopcount<<"\tstacksize="<<xisuStack.size()<<endl;
#ifdef USING_LIST_AS_STACK
		if(loopcount < 10000)
		{
		cout<<"xisu stack=[";
		int showcount = 0;
		for(list<double>::iterator itr=xisuStack.begin(); itr != xisuStack.end(); ++itr)
		{
			++showcount;
			cout<<" "<<(double)(*itr)<<",";
			if(showcount > 20){cout<<"...";	break;}
		}
		cout<<"]"<<endl;
		}
#endif

#ifndef USING_LIST_AS_STACK
		// using stack container
		jieningSquareMatrix tmpsm = smStack.top();
		smStack.pop();
		double tmpxishu = (double)xisuStack.top();
		xisuStack.pop();
#else
		// using list container
		jieningSquareMatrix tmpsm = smStack.back();
		smStack.pop_back();
		double tmpxishu = (double)xisuStack.back();
		xisuStack.pop_back();
#endif
		if(abs(tmpxishu) < jieningMatrix::ZEROVAL)
		{
			cout<<"in doDetComput loop="<<loopcount<<"\t xisu is zero just skip"<<endl;
			cout.flush();
			continue;
		}

		cout<<"in doDetComput loop="<<loopcount<<"\ttop matrix xlen="<<tmpsm.m_xlen<<endl;
		if(3 == tmpsm.m_xlen)
		{
			double tmptmpv = 0.0;
			tmptmpv += tmpsm.m_values[0]*(tmpsm.m_values[4]*tmpsm.m_values[8] - tmpsm.m_values[5]*tmpsm.m_values[7]);	
			tmptmpv -= tmpsm.m_values[1]*(tmpsm.m_values[3]*tmpsm.m_values[8] - tmpsm.m_values[5]*tmpsm.m_values[6]);	
			tmptmpv += tmpsm.m_values[2]*(tmpsm.m_values[3]*tmpsm.m_values[7] - tmpsm.m_values[4]*tmpsm.m_values[6]);	
			tmpv += tmpxishu*tmptmpv;
		}
		else
		{
			int tmpfuhao=1;
			for(unsigned int i=0; i< tmpsm.m_xlen; ++i)
			{	
				if(abs(tmpsm.m_values[i]) > jieningMatrix::ZEROVAL)
				{
					double pushxishu = tmpxishu * ((double)tmpfuhao) * (tmpsm.m_values[i]);
					if(abs(pushxishu) > jieningMatrix::ZEROVAL)
					{
#ifndef USING_LIST_AS_STACK
						// using stack container
						xisuStack.push(pushxishu);
						smStack.push(tmpsm.getSubMatrix(i, 0));
#else
						// using stack container
						xisuStack.push_back(pushxishu);
						smStack.push_back(tmpsm.getSubMatrix(i, 0));
#endif
					}
				}
				tmpfuhao = -tmpfuhao;
			}
		}
		//if(loopcount == 3000) {cout<<endl<<"in doDetComput loop="<<loopcount<<"\tstacksize="<<xisuStack.size()<<endl;break;}
		cout.flush();
	}
	return tmpv;
}

ostream& operator<<(ostream& out, jieningMatrix& mobj)
{
	out<<mobj.toString();
	return out;
}

ostream& operator<<(ostream& out, const jieningMatrix& mobj)
{
	out<<mobj.toString();
	return out;
}

